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Abstract 

A suitable measure for the similarity of shapes represented by param- 
eterized curves or surfaces is the Frechet distance. Whereas efficient algo- 
rithms are known for computing the Frechet distance of polygonal curves, 
the same problem for triangulated surfaces is NP-hard. Furthermore, it 
remained open whether it is computable at all. 

Using a discrete approximation we show that it is upper semi- computable, 
i.e., there is a non-halting Turing machine which produces a monotone 
decreasing sequence of rationals converging to the Frechet distance. It 
follows that the decision problem, whether the Frechet distance of two 
given surfaces lies below a specified value, is recursively enumerable. 

Furthermore, we show that a relaxed version of the Frechet distance, 
the weak Frechet distance can be computed in polynomial time. For this, 
we give a computable characterization of the weak Frechet distance in a 
geometric data structure called the free space diagram. 



1 Introduction 

Suitable distance measures for comparing the similarity of shapes are an impor- 
tant issue in application areas such as computer vision and pattern recognition. 
The distance measure that comes to mind first is the Hausdorff-distance which, 
between two compact sets A, B c M. d , is defined as 

5h(A, B) = max (maxmin lla — 611, maxmin lla — 611) 

v aeA beB beB aeA ' 

where || ■ || denotes the underlying norm on Mr, for example the Euclidean norm. 
So, the Hausdorff-distance considers the maximum distance of a point in one 
set to the other set and it gives reasonable results in many cases. Moreover, it 
can be computed efficiently in two dimensions if A and B consist of sets of line 
segments A 1 5159") or, more generally, in arbitrary but fixed dimension d where 
A and B consist of sets of fc-dimensional simplices |ABG + 03] . 

However, if shapes are modeled by curves or surfaces, there are examples of 
objects having little resemblance but a small Hausdorff distance, see Figure [TJ 
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Figure 1: Nonsimilar polygonal chains with a small Hausdorff-distance 5. 



In these cases the Frechet distance is more appropriate (see [A"G95, 



which is a metric for parametrized geometric objects and was first introduced 
by Frechet for curves [Fre06j and later for surfaces [Fre24j. The idea of the 
Frechet distance is to take into account the "flow" of the curve or surface given 
by its parameterization. A popular illustration of the Frechet distance of two 
curves is the following. Suppose a man is walking his dog on a leash. The man 
is walking on one curve and the dog on the other. Both may stop but not walk 
backwards. Then the Frechet distance is the shortest length of leash allowing 
them to walk on the two curves from start to end. 
Formally the Frechet distance is defined as follows. 

Definition 1. Let /, g be parameterizations of curves or surfaces, i.e., contin- 
uous functions 

f,g: [0,l] fc ^R d , fcG {1,2}, d > k. 
Then their Frechet distance is 



Sf(L9) 



inf 

<r:[0,l] fc ^[0,l] 



sup 

te[o,i] fe 



where the reparameterization a ranges over all orientation preserving homco- 
morphisms. 

In this paper, we assume that the norm || ■ || underlying the definition is 
the Euclidean norm but the results hold for the L\-, and Loo-norm as well, 
and with some modifications can be generalized to other norms, such as any 
Lp-norm, p € N. 

For dimension k = 1 of the parameter space, in particular for polygonal 
curves, the Frechet distance is known to be computable in polynomial time 
AG95]. For two-dimensional surfaces, however, the computation of the Frechet 
distance surprisingly seems to be much harder. In fact, Godau showed in his 
Ph.D. thesis God98] that computing the Frechet distance between triangulated 
surfaces, even in two-dimensional space, is NP-hard. It remained open, however, 
how hard the problem really is, not even its computability could be shown. 

For the special case that both surfaces are simple plane polygons (in a space 
of arbitrary dimension d) Buchin et al. [BBW06J showed that the Frechet dis- 
tance can be computed in polynomial time. 

The first part of this paper contains a partial result concerning the com- 
putability of the Frechet distance. More specifically, we will show that the 
Frechet distance between triangulated surfaces is upper semi- computable, i.e., 
there is a non-halting Turing machine which produces a monotone decreasing 
sequence of rationals converging to the result. It follows that the decision prob- 
lem whether the Frechet distance of two given surfaces lies below some specified 
value is recursively enumerable. 
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The computationally hard part of computing the Frechet distance for di- 
mensions k > 1 seems to be, that according to the definition, the infimum over 
all homcomorphisms of the parameter space has to be taken. For dimension 
one the orientation-preserving homeomorphisms on [0, 1] are "only" the contin- 
uous, onto, monotone increasing functions on [0, 1]. For higher dimensions the 
homeomorphisms can be much "wilder" . 

We tackle this problem in our algorithm by approximating the homeomor- 
phisms by discrete maps which are easier to handle. We do this by first approx- 
imating arbitrary homeomorphisms by piecewise linear homeomorphisms which 
is a known result from topology. These piecewise linear homeomorphisms are 
then approximated by mesh homeomorphisms, i.e., ones that are compatible 
with certain subdivisions of the original triangulations of the parameter spaces. 
Finally, the distance obtained by considering only mesh homeomorphisms can 
be approximated for fine subdivisions by considering the distances at a finite 
number of points. It remains open, whether the Frechet distance between tri- 
angulated surfaces is a computable function in the strong sense. 

The second part of this paper considers a relaxed version of the Frechet 
distance which we will call the weak Frechet distance denoted by S w f- It is 
defined by modifying Definition[T]so that the the reparameterizationer of g which 
is an orientation preserving homeomorphism is replaced by reparameterizations 
a, t of both surfaces which are surjective continuous maps. In the man-dog 
illustration for curves this means that both can choose arbitrary starting and 
ending points, can move forward and backward, and each one has to traverse 
his curve completely. 

We will show that the weak Frechet distance between two given triangulated 
surfaces can be computed in polynomial time. This is done by first considering 
the corresponding decision problem, i.e., the question, whether for given f,g, 
and e > the distance 5 w r(f, 9) < £■ The decision problem can be characterized 
by geometric properties of the so called free space diagram, a data structure 
originally introduced to solve the decision problem for polygonal curves QYG95J. 
This characterization leads to a polynomial time algorithm for the decision 
problem. Since it can be shown that the exact distance must be one of a finite 
set of critical values, we can use the decision routine in a search strategy to 
actually compute the weak Frechet distance. 

2 Semi-computability of the Frechet distance 

We assume that the input to our algorithm are two triangulated surfaces in 
space R d ,d > 2, which are represented by piecewise linear parameterizations 
/, g : [0, l] 2 — > K d . For simplicity, we will denote the surfaces themselves by / 
and g, as well. 

Piecewise linear means that the parameter spaces of / and g are triangulated 
and on each triangle A = (u, v,w), f and g, respectively, are linear maps in the 
sense that f{X\u + A 2 v + X3U1) = Ai/(u) +A2/(«) + Xsf(w) for all X%, A2, A3 > 
with Ai + A2 + A3 = 1 and g has an analogous property. 

We denote the triangulated parameter spaces of / and g by K and L. The 
vertices of the individual triangles have rational coordinates, and the coefficients 
describing the linear maps are rational as well. Thus, a problem instance has a 
canonical finite representation which can be given as input to a Turing machine 
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(we will describe the algorithm in some high level language, however.) 

We will show that the Frechet distance between triangulated surfaces is 
computable in a weak sense according to the following definition which has 
been considered in the research community concerned with the computability 
and complexity of real functions (see, e.g., [WZOO] ). 

Definition 2. A function ip : N — > M is called upper (lower) semi- computable 
iff there is a Turing machine which on input x outputs an infinite, monotone 
decreasing (increasing) sequence of rational numbers converging to <p(x). 

Now we can formulate the first main result of this article: 

Theorem 1. The Frechet distance between two triangulated surfaces in space 
K d , d > 2, is upper semi-computable. 

Theorem [T] immediately implies the following corollary, where (f,g,a) de- 
notes some standard encoding of a triple consisting of two triangulated surfaces 
/ and g, and some rational a > 0. 

Corollary 1. The set {(f,g,a}\ £>f(/,s) < a}, i.e., the decision problem for 
the Frechet distance is recursively enumerable. 

Proof. Consider the Turing machine producing a monotone decreasing sequence 
converging to £>f(/, <?) which exists by TheoremQ] Stop this Turing machine as 
soon as it produces a value less than a. Thus, the algorithm will eventually halt 
for all triples (/, g, a) in the language and it will run forever for the ones not in 
the language. □ 

The computability of Sf in the strong sense of computability theory of real 
functions (see, e.g., [WeiOOj ) remains open, since the sequence produced by the 
algorithm proving Theorem [T] is not shown to effectively converge to (5p (/,</). 
That is, we cannot compute an upper bound on the distance of the value pro- 
duced by the algorithm after k steps to the real value 5f(/, #)■ 

Also observe that Corollary [T] cannot be deduced from Theorem [T] any more, 
if we replace the <-sign in the definition of the set by a <-sign. 

Our proof can be modified to show a weaker form of Theorem [T] for more 
general surfaces. More precisely, if we just assume that the parameterizations / 
and g are computable real functions, it is still correct that there is an algorithm 
producing on input /, g (represented, say, by the Turing machines computing / 
and g) an infinite sequence of rational numbers converging to <5f(/, <?)■ However, 
this sequence is not necessarily monotone decreasing, i.e., the corollary cannot 
be deduced any more. 

2.1 Approximating the homeomorphisms 

We will prove the semi-computability of the Frechet distance by showing that 
homeomorphisms can be approximated arbitrarily closely by so-called mesh 
homeomorphisms, which are compatible with finite subdivisions of the param- 
eter space. 

First, let us recall some standard definitions and notations from topology. 
For a triangulation T let T m denote its m-th barycentric subdivision, where in 
one subdivision step each triangle is subdivided into 6 triangles by the bisectors 
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of its sides. By |T| we denote the underlying space of T, i.e., the set of all points 
lying on triangles of T. In our case \K\ = \L\ = [0, l] 2 . And by mesh(T) we 
denote the maximal diameter of triangles in T which gives a measure for the 
"fineness" of T. 

Now we define mesh homeomorphisms. 

Definition 3. Given two triangulations K and L, a piecewise linear homeomor- 
phism h : \K m \ — > \L n \ is called a mesh homeomorphism if it maps the edges of 
K' m to edge chains of L n , i.e., polygonal chains consisting of edges of L n . 

We will show that any homeomorphism can be approximated arbitrarily 
closely by a mesh homeomorphism. In fact, we need only a weak form of close- 
ness which is defined as follows. 

Definition 4. Let K be a triangulation of the unit square and let h, h' : 

[0, l] 2 — > [0, l] 2 be homeomorphisms on the unit square. Then we define 

d K (h, h!) = max S H (h(A), h'(A)) 

where A E K ranges over all triangles in K and 8h denotes the Hausdorff 
distance. 

Lemma 1. Let K and L be triangulations, a : \K\ — *■ \L\ a homeomorphism, 
m G N, and e > 0. Then there exist n G N and a mesh homeomorphism 
h : \K m \ -> \L n \ such that d Km (a,h) < e. 

Proof. By a theorem from topology (see e.g. chapter 6 in the book by Moise 
}Moi77] ) . a can be approximated arbitrarily closely by a piecewise linear home- 
omorphism, i.e., for all e± > there exists a piecewise linear homeomorphism 
h! : | if | — > \L\ with dK™{<J,h') < e\. We use this fact as a first step, be- 
cause piecewise linear homeomorphisms are much easier to handle than arbitrary 
homeomorphisms. 

We will show that any piecewise linear homeomorphism h' can be approxi- 
mated to any £2 > (in the sense of Definition by a mesh homeomorphism 
h, i.e., dK™{h, h!) < £2- Choosing E\ and £2 so that £1 +£2 < £ then proves the 
lemma, see Figure O 

In order to approximate h! we first show how to find edge chains in L n , for 
some large enough n G N, that are close to the polygonal chains which are the 
images of edges of K m under h! . Then we explain how this can be extended to 
a piecewise linear homeomorphism on the whole parameter space |if m |. 

In fact, the piecewise linear images of the edges of K m under h! form a graph 
G' isomorphic to K m embedded in \L\ where the edges of G' are polygonal chains 
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(a) Construction in the neighbor- 
hood of a vertex 



(b) Edge chain approximating 

7T 



Figure 3: Approximating a piecewise linear homeomorphism by a mesh home- 
omorphism. 

in \L\. We want to modify G' to obtain an isomorphic graph G embedded in 
\L n \ with edge chains of L n as edges that have a distance smaller than £2 to the 
corresponding edges of G', 
We do this in two steps: 

Step 1. We map the nodes of the graph G' and short initial segments of their 
incident edges to nearby vertices of L n and (some of) their incident edges, for 
some suitable n' € N. 

More precisely, for mapping the nodes, we put a small circle of radius r 
around each node h'(v) where r < £2/2 and less than the smallest distance 
between two nodes h'(vi) and h'{v2)- r < £2/2 yields disks of diameter less 
than £2 in which we may move nodes and edges freely, r less than the smallest 
distance between two nodes ensures that disks around different nodes do not 
touch. 

We choose n! S N such that mesh(L n ) < r/2 and for all nodes h'(v) there is 
a vertex w S L n of distance less than r/2 whose degree (in L n ) is greater than 
the degree of h'(v) (in G'). This is possible because mesh(L n ) tends to zero 
for large n' and in the barycentric subdivision the degrees of vertices double 
in each subdivision step after their introduction. For each node v € K m we 
choose h{v) — w for such a w, i.e., satisfying \\h'(v) — h(v)\\ < r/2 < £ 2 and 
deg(h(v)) > deg(v). By this construction each h(v) and all its incident edges of 
L n lie in the disk of radius r around h'(v), see Figure [3](a). 

We start mapping the edges of G' by first choosing edges incident to h(v) in 
L" as initial segments. We do so maintaining the order given by G' , i.e., if the 
edges ei, efc leave the vertex h'(v) of G' in clockwise order, we choose corre- 
sponding edges fi,—,fk 01 L n leaving h(v) in the same order (see Figure E]a). 

Next we cut the edges e\, e/. at the points p±, ■■■,Pk where they first leave 
the disk. Within the disk we connect the free endpoint of each fi with the point 
Pi, i — 1, k by nonintersecting polygonal chains, which replace the original 
polygonal chains from h'(v) to pi. This is possible because we chose the edges 
fi, fk in the same order as the cutting points and within a disk we are allowed 
to move freely, without violating the distance bound £2 (see Figure [3]a). Thus, 
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we replaced all vertices h'(v) by closeby vertices h(v) lying on the mesh L n . 

Step 2. Next, we delete the nodes h(v) together with the incident edges 
from the scenery, leaving a finite set of pairwise disjoint polygonal 
chains which start and end in mesh points of L n . We show that they can be 
^-approximated by edge chains of L n , for suitable n> n'. 

To achieve this, let r\ be the minimum distance between a vertex of a curve 
and a nonincident edge of a curve (possibly the same). For edges sharing an 
incident vertex let 9 be the minimum distance between the intersection points 
of the edges with a circle of radius £2 /2 around their common vertex. 

Then we subdivide L n to L n for a suitable n > n' such that mesh(L n ) < 
mm(i]/2,8/2). Thus, L n is so fine that each polygonal chain it traverses a 
sequence of triangles T v so that T w and even its neighboring triangles are not 
intersected by other polygonal chains (see Figure [21(b)). 

Now, we can approximate each polygonal chain tt by any connected, not self- 
intersecting edge chain a T of L n that lies within T v and connects the endpoints 
of 7r so that the Hausdorff distance between tt and a n is less than e 2 - 

Thus, we showed the existence of a polygonal chain h(e) in L n for each edge 
e of K m which is arbitrarily close to c(e). The chains hie) and the vertices h(v) 
form an embedded graph G isomorphic to G' and, therefore, to K m . h can be 
extended to the interior points of each edge e to form a homeomorphism on e in 
a straightforward manner. Furthermore, the faces of G induce a partition of the 
set of triangles of L n which is isomorphic to the triangulation K m . To extend 
h to a piecewise linear homeomorphism on |if m |, we subdivide each triangle A 
of K m according to the triangulation of the associated set of triangles in the 
partition of L n and extend h to the interior of A correspondingly. □ 

The fact, that homeomorphisms can be approximated arbitrarily closely by 
mesh homeomorphisms implies that the Frechet distance can be approximated 
arbitrarily closely by considering only mesh homeomorphisms and distances 
between the vertices of the corresponding triangulations. More precisely, we 
obtain: 

Lemma 2. 

M/,5) = m f . . , max max \\f(v) - g{w)\\ 

m,n£N h:\K m \— >\L n \ AeKjp vEVa 

where h ranges over all orientation preserving mesh homeomorphisms, K™ is 
the set of triangles in K m , Va are the vertices of A, and MJf,^ is the set of 
vertices of L n that lie in h(A). 

Proof. In order to prove Lemma [5] we define the right hand side of the equation 
in the lemma as the discrete Frechet distance 5dF{f,g) of surfaces and we show 
that it is equal in value to the Frechet distance. 

We first show that the Frechet distance is not larger than the discrete Frechet 
distance. 



Claim 1. 8 F < 6, 



dF- 



7 



Since any mesh homeomorphism is, in particular, a homcomorphism it suf- 
fices to show that for a mesh homeomorphism h : \K m \ — > \L n \ we can bound 
the pointwisc maximum by the maximum taken at vertices, i.e., show that 

max ||/(f)-s(M*))ll < max max \\f(v) - g(w)\\. (1) 

te[0,l] 2 AeK™ vEVa 

To see this, let t € [0, l] 2 be arbitrary, t lies in a triangle A = (vi,V2, 1*3) of 
K m and h(t) lies in a triangle A' = (wx,W2,Ws) of h(A) C L n . Since / and g 
are piecewise linear and K m and L n are refinements of the underlying triangu- 
lations of the parameter spaces, /(A) and g(A') are triangles, as well, namely 
(f(vi),f(v 2 ),f(v 3 )) and (g(wi),g(w 2 ),g(w 3 )), respectively. Consequently, since 
the maximum distance between points of two triangles in 3-space is attained be- 
tween two corners, we have that — g(h(t))\\ < — g(wj)\\ for some 
i, j with 1 < i,j 1 < 3. Taking the maximum on both sides yields equation ([T]). 

Now we show that also the discrete Frechet distance is not larger than the 
Frechet distance. 



Claim 2. For all e > 0, SdF < S F +e. 

The idea is that for any homeomorphism there is a mesh homeomorphism 
arbitrarily close and for the mesh homeomorphism the distance computation at 
vertices comes arbitrarily close to the distance computation on all parameter 
values by sufficient subdivision of the domain complex. 

Let a be a homeomorphism close to realizing Sf, i.e., max ||/(t) — <?(cr(t))|| < 

Sp + ei for some small E\ > 0. By Lemma 1, for any £2 > and any m € N 
there is a mesh homeomorphism h : \K m \ — * \L n \ such that d,K™(cr, h) < £2 ■ 

Let A be some triangle in \K m \ and v one of its vertices. Since d,K m (c, h) < 
£2, for any w £ h(A) C L n there is an x € c(A) with ||ui — x\\ < £2. Using 
t = <t~ 1 (x) and the Lipschitz-continuity of g we get ||g(w) — g(a(t))\\ < c g ■ £2 
for some i€ A where c g denotes the Lipschitz constant of g. 

t and v lie in the same triangle A S K m , so \\v — t\\ < mesh(i^ m ) and, thus, 
\\f{v) — f(t)\\ < Cf ■ mesh(K m ) where c/ is the Lipschitz constant for /. 

Putting everything together and applying the triangle inequality repeatedly 
we obtain 



6 dF < max max \\f(v) - g(w)\\ 

max max II f (v) — q(x)\\ + c a ■ £2 

xecr(A) 

max max || f (t) — q((r(t))\\ + c ■ £2 + Cf ■ mesh(K rn ) 
< S F + £1 + c g ■ £ 2 + cj ■ mesh(K m ). 



Since £i,£2, and, by choosing m large enough, mesh(X m ) can be made arbi- 
trarily small, Claim 2 is shown. 

Claim 1 and 2 yield Sf = SdF which proves Lemma [2] □ 
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2.2 Semi-computing the Frechet distance 

Using Lemma[S]we can now give an algorithm showing the upper semi-computability 
of the Frechet distance between surfaces as claimed in Theorem [TJ This algo- 
rithm will, on input /, g run forever and produce a monotone decreasing se- 
quence of rational numbers converging to <5f(/, g)- 

Algorithm 1: SemiComputeFrechet(/, g) 

Input: Triangulated surfaces f,g, including triangulations K, L of the 

parameter spaces, in a finite description 
Output: A monotone decreasing sequence of rational numbers 
converging to S F (f,g) 

1 set D = oo 

2 forall (n, m) € N X N do 



3 
4 
5 
6 
7 

8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 



generate the barycentric subdivisions K m of K and L n of L 
let E = {ei, efe} be the set of edges in K m 
forall k-tuples {tt±, ...,7Tfc) of simple polygonal chains in L n 
assign to the edge the polygonal chain m for i = 1, ... 
if this assignment results in an orientation preserving 
homeomorphic image of K m then 
set M = 

forall triangles A of K m do 

let Ha C \L n \ be the region in L n assigned to A 
forall vertices v of A and vertices w of Ha do 
| set M = max(M, ||/(u) - g{w)\\) 
end 
end 

set D = min(D, M) 
output D 
end 
end 



do 

k 



19 end 



Line 2 can be realized by some standard enumeration method for pairs of 
integers. The number of fc-tuples of polygonal chains of L n checked in line 5 
is finite. In fact, it is bounded by (l\) k where I is the number of edges in L n , 
which itself is exponential in n, whereas k is exponential in m. But efficiency is 
not the issue here. 

In line 12 we assume that the norm || • || underlying the Frechet distance 
can be evaluated by rational operations. This is correct for, e.g., the L\- or 
Loo-metric but not directly for the Euclidean metric £2- In that case, one 
should rather operate with the square of the distance in line 12 and output 
some suitable rational approximation of y/~D (which is possible) in line 16. 

In line 7 we check whether an assignment of edges in K m to polygonal 
chains in L n results in an orientation-preserving homeomorphic image of K m 
by checking the following three conditions 



the edges on the boundary of | K T 
preserving the orientation 



are mapped onto the boundary of \L r ' 



if a set of edges in K m share an endpoint, the corresponding chains in L n 
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do, as well, 

- other than that, there are no intersection points between two chains. 

Note that checking that the boundary of \K m \ is mapped orientation preserv- 
ingly onto the boundary of \L n \ entails that the mesh homeomorphism is orien- 
tation preserving also on the interior. 

For each pair (m, n) € N x N all mesh homeomorphisms h : K m — ► L n are 
evaluated by Algorithm [U i.e., 

5h,m,n= max max \\f(v)-g(w)\\ 

AG-R"™ uGVa 

(see Lemma [5]) is computed^. 

To see that Algorithm[T]produces values arbitrarily close to Sp(f,g), observe 
that any neighborhood of 5f(/, <?) must, by Lemma [21 contain some value of 
the form Sh, m ,n- The algorithm will eventually encounter that pair (m,n) and 
the subdivision corresponding to h and output 6h, m ,n- By line 15 the output 
sequence is monotone decreasing. Since for all triples (h,m, n), by Lemma 
5h,m,n > £f(/> 9)i l me 15 is justified. Therefore, Algorithm [T] arbitrarily closely 
approximates Sf(J, g) which proves Theorem [TJ 



3 Computability of the Weak Frechet Distance 

In this section we give a polynomial time algorithm for computing the weak 
Frechet distance between triangulated surfaces. The weak Frechet distance is 
a relaxed version of the Frechet distance which does not require the reparam- 
eterizations of the curves or surfaces to be injective. Instead it uses surjective 
continuous maps as reparameterizations. 

Definition 5. The weak Frechet distance between curves or surfaces given by 
continuous parameterizations /, (7 : [0, l] fc — > M. d with k € {1, 2}, k < d is 

<W(/,flO= , mf r it max ||/(a(x)) - g{0(x))\\ , 

a, (3: [0,l] fe ^[0,l] fc x£[0,l]* 

where a and /3 range over all surjective continuous maps. 



As was mentioned before in the man-dog illustration for curves this means 
that both are allowed to choose their start and endpoints, can move forward 
and backward, and eventually must have traversed both curves completely. The 
weak Frechet distance is always less than or equal to the Frechet distance and 
greater or equal the Hausdorff distance: 

S H (f,g) < 6 wF (f,9) < S F {f,g). 

Figure [31 e.g., shows two curves with a small weak Frechet distance and a 
large Frechet distance. For curves, this relation between the distance measures 

1 A more detailed analysis shows that, in fact, it suffices to consider only the pairs 
(m, 2m), m e N 
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Figure 4: Curves with small small weak Frechet distance and their free space 
diagram for parameter e 



is reflected in the free space diagram which is a geometric data structure for 
computing the Frechet distance of polygonal curves [AG95 . 

For the sake of illustration we will shortly repeat these ideas and concepts 
for curves before we extend them to surfaces. For two parameterized curves /, g 
and parameter e the free space is defined as F e (f,g) = {(x,y) | ||/(x) — g{y)\\ < 
e\. The free space diagram illustrates the free space, see Figure |4j where the 
direction from left to right corresponds to the parameterization of / and the 
one from bottom to top to the parameterization of g. The white area inside the 
diagram corresponds to the free space for the indicated value of e. For polygonal 
curves /, g of lengths n and m, respectively, the free space diagram consists of 
nm cells, each one corresponding to the free space between two segments of the 
curves. 

It can easily be seen, that Sp(f,g) < s if and only if there is a monotone 
path in the free space from the lower left corner (0, 0) to the upper right corner 
(1, 1) (corresponding to the reparameterization a in Definition [lj . This fact is 
used in |AG95| to compute the Frechet distance. 

But also the fact that the Hausdorff distance or the weak Frechet distance 
are less than or equal s can be characterized by the free space diagram. Let 
us call the free space or any subset of it extensive if its projections onto both 
parameter spaces are surjective. Then: 

Hausdorff distance less than or equal e is equivalent to the free space being 
extensive. 

Weak Frechet distance less than or equal e means that there is a curve inside 
the free space which is extensive. This fact is used in AG95 to compute a 
restricted version of the weak Frechet distance called the non-monotone Frechet 
distance of polygonal curves. An extensive curve exists exactly if there is an 
extensive connected component of the free space. 

Essentially, this concept and this characterization of the weak Frechet dis- 
tance can be extended to surfaces (yielding a four-dimensional free space dia- 
gram). We will use this idea to find an efficient algorithm for computing the 
weak Frechet distance of triangulated surfaces. 

In fact, the second main result of this article is the following theorem. 

Theorem 2. The weak Frechet distance between triangulated surfaces can be 
computed in 0((m 2 n + mn 2 ) log 2 (mn)) time. 

The algorithms we give will need to compute the intersection of ellipses, 
circles and line segments and compare such intersection points. These intersec- 
tion points can all be described as roots of polynomials of degree up to four, 
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which can be compared and computed exactly using constantly many arithmetic 
operations on rationals, see [ET041 lMPS+06] . 

3.1 Free Space Diagram of Triangulated Surfaces 

As in Section [21 we assume we are given piecewise linear parameterizations 
f,g : [0, l] 2 — * M. d , d > 2, of triangulated surfaces with underlying trian- 
gulations K, L, respectively. All coefficients of the maps and coordinates of 
the triangulations are assumed to be rational. In the following, we will use 
x, y to denote points in parameter space, i.e., x, y G [0, l] 2 . Let F s (f,g) — 
{(x,y) | ||/(x) — ,g(y)|| < e} be the free space diagram of the surfaces. The free 
space diagram lies in the product of the parameter spaces of the surfaces, i.e., 
in the four-dimensional cube. As in the case of curves, we can partition the free 
space into cells. 

Definition 6. Let /, g : [0, l] 2 -> R d , d> 2, be piecewise linear surface param- 
eterizations with underlying triangulations K and L, respectively. A cell of the 
free space is C e (Ax, Al) '■= F e (f,g) n {Ak x Al) for triangles of K and 
Ai of L. A boundary cell is the cell of a triangle and an edge, defined anal- 
ogously. Two cells are called neighboring if they share a nonempty boundary 
cell. 

We will consider the projection of the free space onto the two parameter 
spaces, i.e., the image of F E (f, g) under Proj x : [0, l] 4 — > [0, l] 2 , (x, y) i— > x, and 
Proj L : [0, l] 4 — > [0, l] 2 , (x, y) i— > y. We will compute the projections of the free 
space by computing it for its cells, i.e., we will compute 

Proj^(C £ (A K ,A L )) = {xe A K | Bye A L : ||/(x) - g(y)\\ < e}, 

and analogously for the projection onto the parameter space L. Although the 
free space lies in the product of the parameter spaces, it is defined based on the 
distances in image space. 

The combinatorial structure of the free space is captured by the graph whose 
vertices are the cells of the free space and edges exist between neighboring cells 
which share a non-empty boundary cell. For a free space diagram F e (f,g) we 
denote the corresponding graph by G £ . The graph G e has mn vertices and 
0(mn) edges. For e\ < £2 the graph G ei is a subgraph of G E2 . For varying e, 
G e changes only at finitely many values, namely when a boundary cell becomes 
non-empty. This happens when e equals the minimum distance between points 
on some edge ex of K and on some triangle of L or vice versa: 

C e (e K ,A L )^0 ^ min{\\x-y\\\xef(e K ),yeg(A L )}<E. (2) 

The main idea of our algorithm is contained in the following lemma. 

Lemma 3. The weak Frechet distance between two triangulated surfaces is less 
than or equal e if and only if there is an extensive connected component A in 
the free space for the parameter e, i.e., Projx(A) = Proji,(A) = [0, l] 2 . 

Proof. Let a, (3 be two reparameterizations of / and g, respectively. We define 
a set M a ^ c [0, l] 4 by 

M a , p = {(a(x),/3(x)) |xe [0,1] 2 }. 
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The set M„^ is connected because a and f3 are continuous. M a ^p is extensive 
because a and (3 are surjective. By definition of M a> p the following equivalence 
holds: 

^CF E (/ )5 )« max ||/(a(x))-. 9 (/3(x))|| < e. (3) 

x6[0,l] 2 

Now we show the two directions of the lemma. 

=X If a,/? are two reparameterizations for / and g, respectively, realizing 
weak Frechet distance < e, i.e., max xC [ .i]2 ||/(a(x)) — g(/3(x))|| < e then M a ^ 
is an extensive, connected subset of F £ (f,g), so the connected component con- 
taining M Q)i g is an extensive connected component in F £ (f,g). 

A weak Frechet distance of at most e, however, does not imply that two real- 
izing reparameterizations a, (3 exist, as it is defined as ini a ^ max( x j ||/(a(x)) — 
g(/3(x))|| < e. This can be reformulated as 

Vi?>e 3a, 0: max ||/(a(x)) - g(/3(x))\\ < ■&. 
xe[o,i] 2 

By equation [3] this implies that F#(f,g) contains for all 1? > e an extensive 
connected component. We will show that this property holds also for F £ (f,g). 

As discussed in Section 13.11 the combinatorial structure of the connected 
components of the free space changes only at finitely many values for e, i.e., 
when a boundary cell becomes non-empty. Let 

r\ = argmin{$ | d > e and a boundary cell becomes non-empty in F$(f,g)}. 

Then the combinatorial structure of F# is the same for all $ € [e, rj). 

Let (e n )nGN be a sequence in [e, rf) converging to e. F Cn contains for all n an 
extensive connected component. A free space diagram has only finitely many 
connected components, and therefore there must be at least one extensive con- 
nected component in infinitely many of the F Bn . We claim that this connected 
component is also extensive in F £ . Assume this is not the case, i.e., there is a 
point in parameter space that does not lie in the projection of the connected 
component in F £ . Then because the projection of the free space is a closed set 
this must have already been the case for a small neighborhood [e,rj') C [s,rj), 
i.e., for all but finitely many of the e n ,n G N. This is a contradiction to the 
assumption. 

<=: Let A C F E (f, g) be an extensive connected component of the free space. 
We will construct reparameterizations a and f3 of / and g, respectively, such 
that M atP C A and Proj K (M a ^) = Proj x (A) = [0,1] 2 and Proj L (M Qj/3 ) = 
PiO] L (A) — [0, l] 2 . This will imply that a and j3 realize a weak Frechet distance 
less than e. 

We do this in two steps: first we define for each cell C £ (Ak, A^) contained in 
the connected component A two parameterized surface patches Mk(Ak,Al) 
and M L (A K ,A L ), which fulfill Proj K (C E (A*, A L )) = Proj K (M K (A K , A L )) 
and Proj L (C £ (Aif , A^)) = Proj £ (Afi(A^, Al)). By a surface patch we mean 
a two-dimensional manifold with boundary in R 4 . Then we show how to "glue" 
two such parameterized surfaces patches together. Glueing together the surface 
patches of all cells contained in the connected component will give the desired 
reparameterizations. 
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Defining Parameterized Surface Patches in each Cell For each cell 
C s (Ak, Al) of the free space we define two surface patches: 

M K (A K ,A L ) = {(x )fl - 1 (n(/(x), 5 (A i ))))|xeProj A: (C e (A J f,A i ))} ) 
M L (A K ,A L ) = {(.r 1 (n(g(y),f(A K )),y))\yePToi L (C £ (A K ,A L ))}, 

where 

n(x, A) = argmin \\x — y\\ 

is the nearest neighbor of a point x in the set A. 

These surface patches have the following properties: 

1. M K (A K ,A L ) U M L {A K ,A L ) c C e {A K ,A L ). 

2. Vio- ]K {C e {A Kl A L )) = Proj K (M K {A K ,A L )) and 
Proj L (C £ (A K ,A L )) = Pvoj L (M L (A K ,A L )). 

3. M K (A K ,A L ) n M L {A K ,A L ) ± 0. 

4. If the shared boundary cell of the cells C e (Ak, Al) and C e (Ax, A' L ) is 
non-empty, then M L (A K , A L ) n M L (A K , A' L ) ^ 0. 

If the shared boundary cell of the cells C e (Ak, Al) and C S (A' K , Al) is 
non-empty then M K (A K , A L ) n M K (A' K , A L ) ^ 0. 

5. Parameterizations olktPk and ctLiflh exist which fulfill Mk(Ak, Al) = 
M aKt0K and M L {A K , A L ) = M aLi0L . 

For the first property consider a tuple (x, y) € Mk(Ak, Al). By definition 
the point x is in Pto] k {C e {Ak , A^)), i.e., there is a point y' s.t. (x, y') e 
C £ (Ak, Al). But then also (x, y) e C s (Ak, Al) because y was chosen based 
on the nearest neighbor in image space. 

The second property holds by definition, that is, Mk{Ak 1 Al) is defined by 
choosing a tuple for each point x € Proj x (C e (A/f, A^)). 

The third property is equivalent to the existence of two points in Ak and 
Al whose images under / and g, respectively, are each others nearest neighbor, 
i.e., a tuple (x, y) e Aif x A^ such that g(y) = n(/(x), <?(Al)) and /(x) = 
n (#(y): /(Ak)). Points fulfilling this condition are points whose images have 
minimal distance in the free space cell, i.e., a tuple (x, y) such that ||/(x) — 
ff(y)|| < ||/(x') - g(y')\\ for all tuples (x',y') G Ak x A l . Such a tuple exists 
because Ak x Al is compact. 

The fourth property holds because Ml(Ak,Al) and Ml(Ak, A' L ) coin- 
cide on the intersection of C £ {Ak 1 Al) and C e (Ak 7 A' L ), and analogously for 
M K (A K , A l ) and M K (A' K , A L ). 

The fifth property is fulfilled by the following natural parameterizations. 
Let ok and ol be homeomorphisms from [0, l] 2 to Ak and A^, respectively. 
Then a K = o K and (3 K = .9 _1 ° n(-,g(A L )) o f o a K and a L = f^ 1 o 
n(-,/(Ajf)) o g o ol and (3l — &l are parameterizations of Mk{Ak 7 Al) and 
Ml{Ak, Al), respectively. These reparameterizations are continuous because 
the identity map, /, g and the nearest neighbor function on convex sets in R d 
are continuous. 
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Figure 5: Glueing together two parameterizations 



Now we define Ma to be the union of all surface patches of cells in the 
connected component A, i.e., 

M A = |J (M K (A K ,A L )UM L (A K ,A L )). 

C e (A K ,A L )cA 

The first two properties of the surface patches Mk(Ak , A/,) and Ml{Ak, Al) 
imply the following properties of Ma' 

1. M A c A 

2. Proj^(A) = Pto) k (M a ) and Proj L (A) = Proj L (M A ). 



Glueing together the Parameterizations of the Surface Patches We 

first show that two reparameterizations can be joined to one and then obtain 
parameterizations a a, Pa for Ma by successively glueing together parameteri- 
zations of the surface patches in Ma- 



Claim 3. Given parameterizations a, f3, a',/3' with M a $ n M a i ^ 0. Then 
there exist parameterizations a" , (3" such that M at p U M a i^i = M a >t t 0". 

Proof of Claim 3: Let p = (x p ,y p ) be an intersection point of M Q|j g and 
Ma^pi. Let x ,xi be such that p = (a(x ), /9(x )) and p = (a'(xi), /3'(xi)). 

For glueing together the parameterizations a,j3,a',/3', we partition the pa- 
rameter space as illustrated in Figure [5] (a). Choose a point c in parameter 
space and radii r, r' such that r < r' and the disk of radius r' around the point 
c is contained in [0, l] 2 . Let C, C be the circles of radius r,r' around c and 
D, D' the (closed) disks of radius r, r' around c. Let R be the half-open ring 
with radii r, r' around c and let S = [0, 1] 2 \D'. 

We will glue the parameterizations a, (3 and a',/3' along the circle C On 
the circle a" , (3" will be constant, outside they will be defined based on a',/3' 
and inside based on a, (3. 

We define homeomorphisms on the parts of the partition as illustrated in 
Figure [5] (b). Let o\ be a homeomorphism from D to the unit cube. Let 02 be 
a homeomorphism from R to the unit cube without the point xq. Choose o\ 
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and <72 s.t. they coincide on the circle C . Let 03 be a homeomorphism from S 
to [0, l] 2 \{xi}. Now we define a", (3" : [0, l] 2 -> [0, l] 2 as follows: 

'a(cri(x)), x£fl f/3((Ti(x)), x e D 

„ = . a(cr 2 (x)), x e R „, , = I /3(cr 2 (x)), xefi 

S x p , xeC" P 1 ' |y p , xgC' 

a'(0- 3 (x)), xeS U'Mx)), xeS 

a", (3" are continuous and fulfill M a ,£ U M a > t pi = M a ii t pn by construction. 
They are continuous inside D,R, and D as concatenations of continuous func- 
tions. On the circle C they are continuous because o\ and 02 coincide on C. 
They are continuous on the circle C because 02 (x) tends to Xo for x tending to 
x' G C and therefore a(o2(x)) and /3(o2(x)) tend to x p and y p , respectively, 
for x tending to x' E C . Similarly, 03 (x) tends to xi for x tending to x' e C" 
and therefore a'((T3(x)) and /3'(o3(x)) tend to x p and y p , respectively, for x 
tending to x' £ C'. Thus, Claim 3 is proven. 

We can glue together all surface patches in Ma because they are connected 
by intersection - properties 3 and 4 of the surface patches. By glueing together 
all surface patches in Ma we get reparameterizations a a, Pa for /, g which re- 
alize a weak Frechet distance between / and g which is less than e: aA and 
(3a are continuous because the parameterizations of single surface patches are 
continuous and the glueing of parameterizations preserves continuity, ayi and 
(3a are surjective because the projection of M aA ^ A = Ma equals the projec- 
tion of A - property 2 of Ma - and this was assumed to be [0, l] 2 . Finally, 
||/(o;^(x)) — <7(/?a(x))|| < e holds for all x S [0, l] 2 because M aA ^ A lies com- 
pletely in F e (f,g). 

□ 



3.2 Decision Problem for the Weak Frechet distance 

By Lemma [31 deciding if the weak Frechet distance is less than some given value 
e is equivalent to determining if there is an extensive connected component in 
the free space for the parameter e. We will do this algorithmically in two steps: 
First we compute all connected components of the free space, secondly we test 
whether there is an extensive one. 

For the first step we compute the combinatorial graph of the free space as 
described in the previous section in 0{mn) time. 

In the second step we need to determine if the projection of a connected 
component A completely covers both parameter spaces. This property for, e.g., 
the parameter space K is is equivalent to each triangle A^- in K being covered 
by the projection of A. Equivalently, /(A^) in image space must be completely 
contained in the union of the e-neighbor hoods ff(Ai) © B £ of all triangles 
such that the cell C s (Ajc, Al) lies in A. 

We can reduce this to a two-dimensional problem by intersecting the e— 
neighborhoods g{Ai) © B £ with the plane Ea k which contains the triangle 
/(Ak). I.e., we decide for each triangle A#- in K whether 

/(Ajc)C |J ((g(A L )®B e )nE AK ). (4) 

A L £L 
C E (A K ,A L )eA 



16 



This property can be decided by computing for each triangle Ak the arrange- 
ment of the triangle f(Ax) and the set of boundary curves 
E\ K ) for all triangles g L with C e (Ax, Al) G A where d denotes the 
boundary operator. Then we apply a line sweep algorithm to the part of the 
arrangement contained in the triangle f(Ax)- Property (jl]) holds if and only if 
there is no empty cell in this part of the arrangement. 

For the metrics L\ and L^, the e-neighborhoods can be described by lin- 
ear equations. For the Euclidean metric the e-neighborhoods are described by 
quadratic equations: the e-neighborhood of a triangle in 3-space is the union 
of e-balls around the vertices, cylinders of radius e around the edges and a 
triangular prism of height e. The boundaries of the intersections of these e- 
neighbor hoods with a plane are the union of a constant number of half-ellipses, 
circles or half-circles and straight line segments. The event points of the sweep 
over the arrangement are all intersection points of a boundary curve with the 
triangle f(Ax) and intersection points between two boundary curves that lie 
inside f{A K ). 

The overall time complexity of the second step of the algorithm adds up as 
follows: For each triangle Ak we need to decide by which connected components 
of the free space it is completely covered. For one connected component, the 
arrangement for Ak has size 0{l 2 ) and can be swept in time 0{l 2 \ogl) where 
I is the number of cells projecting onto Ak- Let s be the number of connected 
components and Zj, 1 < i < s, the number of cells projecting onto Ak in the 
ith connected component. Because the connected components of the free space 
are disjoint the sum of the li, 1 < i < s, is at most m. Thus, for each triangle 
Ak € K we can determine in 0(X)i=i( m l l°g TO i)) € 0(m 2 logm) time, by 
which of the connected components of the free space it is covered. We have to 
apply the same procedure with K and L exchanged. This yields a total runtime 
of 0(nm 2 logm + mn 2 logn) for all triangles in both parameter spaces. 

Summarizing, we obtain the following algorithm for the decision problem of 
the weak Frechet distance. 
Algorithm 2: Decide WeakFrechet(/, g, e) 

Input: parametrized triangulated surfaces /, g, e > 
Output: Is 5 W F(f,g) < e? 

1 compute the graph G e of the free space diagram F e (/, g) 

2 forall connected components A of G s do 

3 forall triangles A in either parameter space do 

4 using line sweep decide whether A is completely covered by the 
component A 

end 

6 end 

7 output true if a connected component covering all triangles has been 
found, else output false 

Lemma [3] and the analysis above yield the following theorem: 

Theorem 3. Algorithm® decides whether the weak Frechet distance between 
two triangulated surfaces with n and m triangles, respectively, is less than a 
given parameter e in 0{nm 2 logm + mn 2 logn) time. 
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3.3 Computing the Weak Frechet distance 

The decision algorithm can be extended to a computation algorithm by searching 
a set of critical values for e as is done for the Frechet distance between polygonal 
curves |AG95| . We first characterize these critical values and then show how to 
compute the weak Frechet distance by searching over the set of critical values. 

By Lemma [3] the weak Frechet distance corresponds to an extensive con- 
nected component of the free space. Thus, the weak Frechet distance equals a 
value £ if and only if the free space F e contains such a connected component and 
it does not for any smaller value of e. We can characterize the critical values as 
follows. 

1. The combinatorial structure of the connected components changes 

2. The projection of a connected component covers a parameter space com- 
pletely whereas for smaller values of e this is not the case although there 
is no combinatorial change. 

We now describe this subdivision of the critical values more precisely: 

Type 1 The combinatorial structure of the free space changes when a bound- 
ary cell, i.e., the cell of an edge and a triangle, becomes non-empty. This 
happens for e equal to the distance of the edge and the triangle defining the 
boundary cell, as stated in equation ([2|). There are 0(mn) such values, each of 
which can be computed in constant time. 

Type 2 Apart from the combinatorial changes, the projection of a connected 
component grows monotonously with increasing e. We distinguish the critical 
values of this type by the last points in a certain neighborhood covered by the 
projection, i.e., points x € [0, l] 2 which are in the projection of the free space 
F e (f,g) and for any e' < e there is a neighborhood of x which is disjoint to 

F £ '(f,g)- 

A last point covered may be either a vertex (Type 2a) , an interior point of an 
edge (Type 2b) or an interior point of a triangle (Type 2c) of the triangulation 
of the parameter space, or - if the triangle are not in general position - the last 
points covered may be a segment (Type 2d). 

Type 2 a If the last point covered is vertex, then its image (under / or g) 
lies on the boundary of the e-neighborhood of a triangle. Thus, its image has 
distance e to the triangle. We can compute these critical values by computing 
for all triangles of the one surface their distance in image space to the vertices 
of the other surface. Each such value can be computed in constant time and 
there are 0(mn) such values. 

Type 2 b If the last point covered lies in the interior of an edge of one surface, 
then the boundaries of ^neighborhoods of two triangles of the other surface 
intersect in the image of that point. That is, the point has equal distance to 
two triangles. There are 0{mn 2 +m 2 n) such critical values for both parameter 
spaces. 
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Type 2 c If the last point covered lies in the interior of a triangle, then the 
boundaries of e-neighborhoods of three triangles intersect in the point. That is, 
the point has equal distance to three triangles. There are 0(mn 3 + m 3 ri) such 
critical values for both parameter spaces. 

Type 2d If the triangles are not in general position it may also happen that 
the last points covered form a line segment. A detailed analysis shows that this 
can occur only if two edges, or two triangles, or an edge and a triangle in image 
space are parallel and e is the distance between such a parallel pair. There are 
0(m 2 + n 2 ) such values, each of which can be computed in constant time. 

Summarizing, we have 0(m 2 n + n 2 m) critical values of types 1, 2a, 2b, and 
2d, each of which can be computed in constant time. Furthermore, we have 
0(m 3 n + mn 3 ) critical values of type 2c. 

We could now compute the weak Frechet distance by computing and sorting 
the critical values and then doing a binary search over the set of critical values, 
solving Algorithm [2] in each step. The runtime of this algorithm would be dom- 
inated by the sorting of the critical values which takes 0((m 3 n + mn 3 ) log(mn)) 
time. 

We can improve the runtime applying parametric search [Meg83| to the 
0(m 3 n + mn 3 ) critical values of type 2c where the boundaries of three projected 
cells intersect. For the 0(m 2 n + inn 2 ) other critical values we first run a binary 
search involving the decision algorithm[2]in each step. This part takes 0((m 2 n+ 
mn 2 )log(mn)) time. 

More precisely, we have the following algorithm. 

Algorithm 3: ComputeWeakFrechetDist(/, g) 

Input: Parametrized triangulated surfaces /, g 
Output: S wF (f,g) 

1 Compute the set C\ of critical values of type 1, 2a, 2b, and 2d 

2 Sort the set C\ 

3 Do a binary search over the set C\ applying Algorithm [5] in each step and 
branching to larger values if the answer is false and to smaller values if 
the answer is true 

4 Let Ci be the smallest value in G\ for which Algorithm [5] outputs true 

5 Do a parametric search over the set Ci of critical values of type 2c in the 
interval (cj_i,Cj], again applying Algorithm [5] in each step 

6 Output the smallest value for which Algorithm [2] outputs true 

For the parametric search, instead of giving a parallel algorithm for the 
decision problem, it suffices to use any parallel algorithm whose critical values 
include the critical values of the decision problem. In our case, we can use 
a parallel comparison-based sorting algorithm for sorting lexicographically the 
intersection points of the boundary curves of e-neighborhoods of triangles of the 
one surface intersected with a plane containing a triangle of the other surface. 
Formally, we sort all points of the form 

a((A;ffiB E )n£ A ) n d({A' 2 (BB e )nE A ), 

where A is the image of some triangle in one parameter space and , A' 2 images 
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of triangles of the other parameter space. These are 0(m 2 n + n 2 m) many values 
and their critical values include the critical values of the decision algorithm. 

More specifically, we use Cole's variant of parametric search based on sort- 
ing |Col87| . His technique yields a runtime of 0((k + Tdec) log A;) where Td ec is 
the runtime of the decision algorithm and k is the number of values to be sorted. 
In our case, k S 0(m 2 n + n 2 m) and Xd ec € 0(m 2 n log n + n 2 m log to). Thus we 
get a total runtime of 0((m 2 n + mn 2 ) log 2 (tow)), which proves Theorem [5] 
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